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Abstract 

We use lattice Boltzmann simulations, in conjunction with Ewald summation meth- 
ods, to investigate the role of hydrodynamic interactions in colloidal suspensions of 
dipolar particles, such as ferrofluids. Our work addresses volume fractions (j) of up to 
0.20 and dimensionless dipolar interaction parameters A of up to 8. We compare quan- 
titatively with Brownian dynamics simulations, in which many-body hydrodynamic 
interactions are absent. Monte Carlo data are also used to check the accuracy of static 
properties measured with the lattice Boltzmann technique. At equilibrium, hydro- 
dynamic interactions slow down both the long-time and the short-time decays of the 
intermediate scattering function S{q^t)^ for wavevectors close to the peak of the static 
structure factor S{q), by a factor of roughly two. The long-time slowing is diminished 
at high interaction strengths whereas the short-time slowing (quantified via the hydro- 
dynamic factor H(q)) is less affected by the dipolar interactions, despite their strong 
effect on the pair distribution function arising from cluster formation. Cluster forma- 
tion is also studied in transient data following a quench from A = 0; hydrodynamic 
interactions slow the formation rate, again by a factor of roughly two. 

1 Introduction 

The tendency of ferromagnetic colloidal particles to form aggregated structures, stabilized by 
anisotropic dipole-dipole interactions, was insightfully discussed by Pierre-Gilles de Gennes 
and Philip Pincus in 1970 [T]. This tendency is central to the structure f2], and hence the 
phase equilibria [3] HJ [5] and dynamics [6^ 7J, of magnetic colloids in organic solvents (fer- 
rofluids) . The earliest direct observations of chain-like structures were obtained by Hess and 
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Parker using electron microscopy [8] but remarkably the quantitative experimental study of 
strong pair correlations began only in 2003 with cryogenic transmission electron microscopy 
studies by Philipse and co-workers P [TOl [H] . This delay in confirming the classic predic- 
tions of de Gennes and Pincus partly refiects the extreme opacity of most ferrofiuids, which 
precludes both direct microscopy of bulk phases and light scattering as methods to elucidate 
structure. Alongside X-ray and neutron scattering, these two methods have been central to 
the widespread progress made in understanding other forms of colloidal aggregation since 

1970 [laiiaiii]. 

The same experimental difficulties have also made it hard to study structural relax- 
ation, which in macroscopically isotropic ffuids at equilibrium can be quantified by the time- 
dependent correlator (directly accessible in inelastic scattering experiments, where available) 

[13 [IS] 

1 ^ 

S(^^t) = N ^ exp{zg.[r,(t)-n(0)]} (1) 

j,k=l 

where q = \q\. Difficulties with scattering methods mean that many relaxation studies on 
ferroffuids have been limited to strictly q = properties such as frequency dependent bulk 
magnetic susceptibilities [IZlIIH] or magnetoviscous and rheological properties [19j. 

It is increasingly possible for computer simulation methods to fill the gaps in our exper- 
imental knowledge of structural relaxation in complex fiuids [20l ETJ [22l [23l EH [25] . How- 
ever, for ferrofiuids, such simulations have previously been oversimplified in their neglect of 
hydrodynamic interactions between particles. These interactions are mediated by the inter- 
vening solvent - an essentially incompressible, Newtonian fiuid of viscosity r] and density p. 
Previous simulations using molecular dynamics (MD) [26l [271 ESI [29], Monte Carlo (MC) 
[30l [3ll [32] , and Brownian dynamics (BD) j7l[33l[3l] of dipolar fiuids, while each capable of 
generating the correct equilibrium statistics governed by the Boltzmann distribution, are all 
compromised by either the neglect or incomplete treatment of thermal noise and many-body 
hydrodynamics. In particular BD, while capturing the overdamped, diffusive dynamics of 
an isolated particle (with diffusion constant Dq = kBT/dnrja^ a being the particle radius) 
fails to correlate the Brownian motion of two or more particles in the correct manner. As a 
result, S{q^t) will have the wrong time dependence. 

In this work, we present simulation results for dipolar colloids generated using the lattice 
Boltzmann (LB) method. This approach treats full hydrodynamic interactions, at least for 
the colloids simulated here, which have a sufficiently repulsive soft-core potential to ensure 
that particles do not approach one another too closely. (The presence of such a potential 
avoids large hydrodynamic lubrication forces in thin ffuid films between two solid particles in 
close contact, whose treatment within LB is possible, but costly [35].) Alongside LB, other 
methods that fully treat hydrodynamics include force methods [36], Stokesian dynamics 
(SD) and accelerated Stokesian dynamics (ASD) [371 [38], stochastic rotation dynamics 
[39] . (Dissipative particle dynamics does so also, but without proper control of noise terms 
as required here [40j.) Few of these methods have yet been used to treat systems with 
long-range interactions, although an SD method has been applied to ferrofiuids in shear 
fiow 1^ [42], and ASD was recently used to address charge-stabilized colloids at low ionic 
strength [16]. We are aware of no systematic application of these methods to the equilibrium 
dynamics of dipolar fiuids in three dimensions. 
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This paper is organized as follows. In section [2] we outline the numerical methodology, 
and in section |3] present results for equilibrium structure as calculated by BD, LB, and MC 
methods. (These should be, and nearly are, identical.) Then in section [3] we present results 
for S{q^ t) and its orientational analogs, focusing on a like-for-like comparison between BD 
(no hydrodynamics) and LB (full hydrodynamics). In section [5] we present an analysis of 
transient behavior, addressing the time evolution of the static structure including an analysis 
of cluster statistics. Finally in section [6] we give our conclusions and discuss prospects for 
future work. 

2 Simulation methods 

We use the lattice Boltzmann method for a fluid incorporating spherical solid particles 
[35l|4^. In this method, the density, momentum, and stress in the fluid are associated with 
various moments of a kinetic distribution function /(Q,r), deflned at each site r on a 3D 
lattice and acting on a space of 19 discrete velocities q that connect neighboring sites in one 
timestep (including the null velocity). Setting the lattice parameter, timestep, and mean 
fluid density p = {Y,i f) all to unity deflnes a set of LB units that we use in results quoted 
below. (Of course, many physical quantities can be expressed in dimensionless forms, from 
which these units cancel out.) 

Each of our spherical colloidal particles has a hard-core radius a = 2.3 and resides oflF- 
lattice; its surface then cuts the lattice bonds at a set of links at which fluid and particle 
interact by momentum transfer. This transfer is achieved by a 'bounce-back on links' algo- 
rithm j35] . The force and torque on a colloidal particle are found by summing contributions 
across the boundary links; particle velocities and angular velocities are then updated via a 
standard molecular dynamics routine [351 143]. The thermal noise, responsible for colloidal 
Brownian motion, is generated entirely within the fluid and is fully included in the description 
of the fluid momentum [44] , using a method reported previously [45j . Momentum transport 
then causes the random forces (and torques) felt by diflFerent particles to become correlated, 
in accord with the fluctuation-dissipation theorem which relates random forces to the matrix 
of particle mobilities Mij^fs- This matrix obeys Mija(3 = dva^i/dff3j where Vi is the velocity 
of particle i in response to a force fj on particle and are cartesian indices. In the 
presence of hydrodynamic interactions, Mijaf3 is not diagonal in particle indices, but has 
long-range correlations. In contrast to SD-based methods, LB avoids explicit computation 

of Mijc,f3- 

Each colloidal particle has the same mass m = pvQ as the nominal volume of fluid 
that it displaces; '^o = 47ra^/3 is the volume of a colloidal particle. The fluid viscosity 
is set to T] = 0.025 and the temperature to ksT = 5 x 10~^; these choices represent the 
best compromise we have found between numerical accuracy and efliciency for the systems 
under study here. (For a discussion of parameter optimization in LB see [22l E^.) Based 
on these parameters, the particle velocity relaxation time is = m/Gnija 47, and the 
time scale for fluid momentum to equilibrate around a particle is = o^p/rj 212. The 
single-particle diffusivity is i^o = ksT/GTrrja 4.61 x 10~^, giving a diffusive timescale of 
td = cl^/Dq 1.15 X 10^. This offers a reasonably wide domain in which to study short time 
diffusion, even in the case of strongly interacting particles (as here) where the short-time 
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diffusion time window closes off at t > TDih/aY with h a typical surface-to-surface separation 
between adjacent particles. Our choices for other (reduced) simulation parameters ^ to be 
defined below - correspond to typical values for real ferrofluids. Time-dependent results are 
expressed in units of r^^ which means that for t ^ t^^Tj^^ they will be directly relevant to 
almost any real ferrofluid system. 

To avoid computationally expensive lubrication contacts between particles (as mentioned 
in section [T]) we introduce a short-range, soft-core repulsion U^^ acting at separations beyond 
the hard-core radius a = 2.3. (The latter coincides with the hydrodynamic radius [35l [43].) 
We tested various options for U^^ and found that a satisfactory choice is 

Usc(^^ _ f Uo{h) - Uoihc) -{h- K)dU^/dh\h=h. < h < 
^ ^ [ h> he 

where Uo{h) = kBTa/h^ h = r — 2a is the surface-to-surface separation for spheres whose 
centers are r apart, and he is a short-range cutoff. We set the cutoff separation he = 1.2, 
which in terms of the hard-core diameter is hc/2a 0.26. This comprises a truncated and 
shifted inverse power law potential; note that the interaction force remains divergent at 
contact {h = 0). Although structural and magnetic properties of ferrofluids are known to 
be quite insensitive to the choice of short-ranged repulsive potential [29j, our choice could 
qualitatively describe the effects of a screened electrostatic repulsion. The case of steric 
stabilisation is more complex since the polymer layer will have, in addition, a direct effect 
on the hydrodynamic forces. 

The introduction of U^^ reduces discretization errors in the noise forces which become 
acute when fluid nodes are excluded from the space between particles. (Since absence of fluid 
entails absence of noise, particles that are too close to one another to have fluid nodes between 
them effectively feel a reduced temperature.) In fact this issue, rather than avoidance of 
lubrication contacts per 5e, prevents efficient use of a shorter-range soft-core potential than 
the one we have chosen; the effect of using such a potential is to exaggerate the peak in the 
pair distribution function g{r) for particles at close contact (as though particles become colder 
at close separations). Even with our choice of [/^^, this deviation remains visible in the data 
of section [3] for the largest dipole strength used, but this is considered acceptable. Indeed, in 
the current state of the art for LB, errors of several percent from this and other sources such 
as shape discretization remain unavoidable. To reduce these further is straightforward in 
principle: one simply increases a. However, for a fixed number of particles at a given volume 
fraction 0, the system volume must then be increased as a^, with a further increase in the 
run time to reach t^. Thus a factor 2 increase in a, as would be required to give a worthwhile 
reduction in discretization errors, entails a 32-fold increase in computational resource. 

The total colloid-colloid pair potential in our simulations is the sum of short-range soft- 
core (sc) and long-range dipolar (d) contributions: 

Uij = U'''{hij) + U''{rij,Si,Sj). (2) 

Here, Vij is the center-center separation vector for particles i and j, Si denotes a unit vector 
pointing along the dipole of particle z, and hij = r^j — 2a where r^j = |r^j|. The dipole-dipole 
interaction potential is written 



U%r,j,s,,Sj)=8XkBTa' 



(3) 
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where Vij = Vij/rij is a unit vector pointing along the center-center separation vector. The 
interaction parameter A is a dimensionless, dipolar coupling constant defined such that two 
colloids at hard core contact (r = 2a), with dipoles mutually aligned and parallel to r, 
have = — 2A/c^T. This 'nose-to-tail' parallel conformation is of lowest energy; the other 
minimum-energy conformation is 'side- by-side' antiparallel, for which = — A/c^T (with 
corresponding energy maxima on reversing the direction of one dipole). In experimental 
ferrofluids A values up to :^ 4 are readily available; much larger ones can be achieved but 
are relatively unusual P, [101 E] • 

All of our simulations were performed in a cubic box of side L with periodic boundary 
conditions for both fiuid and particles. A standard Ewald summation technique was de- 
ployed to handle the long-range aspect of the dipolar interactions ^ . The Ewald sum 
computes dipolar forces directly in real space for particle pairs with separation r^j < rc, 
a cutoff distance, and deals with the remainder in reciprocal space. To allow a parallel 
implementation using domain decomposition, which is important given the relatively high 
computational requirements associated with the LB fiuid [43j, we chose Tc = 16 in all cases. 
This is small enough that some parallelism is possible, but not so small that the number 
of terms required in the reciprocal space sum for acceptable accuracy becomes unwieldy. A 
convergence factor [47] of k = 5/2rc = 0.15625 was chosen, with wavevectors k = (27r/L)n 
with |n| < 8 and 16 for L = 64 and 128, respectively. These parameters were optimised 
using established methods [18]. Finally, following normal practice, the Ewald sum boundary 
condition at infinity was chosen to be "conducting", i.e., the dipolar (magnetic) susceptibil- 
ity of the surroundings is infinite. This removes a zero-wavevector, bulk-magnetization term 
from the Hamiltonian which, in the opposite extreme of vacuum boundary conditions, can 
lead to unphysically small polarised domains and hence slow down simulation convergence 

m- 

For comparison with our LB results, a BD algorithm was set up within the colloidal MD 
module of our LB code, deploying standard techniques to generate the required Langevin 
dynamics with independent noise acting directly on each particle. The inertia of the par- 
ticles is retained but the many-body hydrodynamics are replaced by a Stokes drag that is 
independent of the location of other particles. By setting exactly the same value for r^^ 
we thus create an algorithm which diflFers from LB solely in the omission of many-body 
hydrodynamics . 

To further validate both codes, and to monitor the achievement of Boltzmann equilib- 
rium, canonical Monte Carlo simulations were performed in a cubic simulation cell with 
periodic boundary conditions applied [47j. The long-range dipolar interactions were han- 
dled using the Ewald summation with conducting boundary conditions, a convergence factor 
K.L = 5.6, and wavevectors k = {2Ti/L)n with \n\ < 6. The maximum translational and 
orient ational displacement parameters were adjusted independently to give acceptance rates 
of approximately 20% and 50%, respectively; it is efiicient to employ low acceptance rates 
for translations of particles with hard cores, due to the possibility of rapidly identifying 
overlaps. For each state point considered, we performed equilibration runs of 2 x 10^ MC 
cycles, where one MC cycle consisted of, on average, one attempted translation or rotation 
per particle. Production runs consisted of 5 x 10^ MC cycles. 

Our LB work was performed primarily on lattices of size V = 64^ or V = 128^. The 
colloid volume fraction is given by = Nvq/V^ where is the number of colloids. For the 
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larger system size at a volume fraction (f) = 0.10, a run of ^ 10^ timesteps required 56 hours 
on 64 cores of a cluster of 3GHz Intel dual-core processors. For resource reasons most of our 
results on fully equilibrated samples {t > 2-3 x 10^ "ibTu at A = 8) concern the smaller V. 
BD and MC runtime requirements were modest in comparison. 

3 Static structure 

In this section we confirm that our LB algorithm generates, to acceptable accuracy, the 
Boltzmann distribution for thermal equilibrium properties, as does our BD code, and that 
both are in agreement with MC data. This is of course necessary if the dynamical data in 
subsequent sections are to be trusted. Such agreement is not automatic but is in fact a very 
demanding test of the ability of our LB algorithm to generate correlated noise forces satisfying 
the fiuctuation dissipation condition. In particular, without adopting the methods of [45] 
(in which noise terms are applied to not only the hydrodynamic but also the local modes of 
the fiuid degrees of freedom) we would not be confident of achieving such agreement. Even 
with these methods, LB parameter values must be carefully chosen to maintain acceptable 
performance. We have already mentioned that correct treatment of noise was the limiting 
factor in our choice of U^^] it also prevents us using a much larger ksT value (which sets the 
intrinsic noise level in the LB fluid) and/or a much smaller viscosity rj. Either step could in 
principle drastically reduce the run time required to reach the basic timescale td for colloidal 
diffusion. Thus the control of discretization errors within the noise sector currently remains 
the primary efficiency bottleneck in the application of large-scale LB simulations to colloidal 
diffusion. 

3.1 Energy equilibration 

Table [T] shows time- averaged energy data for various simulation runs. Those for V = 64^ have 
fully equilibrated (run times > 30td) as judged by convergence of the energy parameters to 
their long-term averages. The results indicate excellent agreement between LB and the other 
methods at A = and 4, and adequate agreement (errors of less than 5%) at A = 8. Some 
runs with the larger system size, V = 128^, are also reported in the Table. These have run 
times ^ 7td and while their energies appear to have saturated, those at A = 8 are showing 
continued structural evolution by other measures (such as cluster statistics; see section [5]). 
Accordingly, the reported energy discrepancies for these runs may include systematic errors 
arising from incomplete structural equilibration, and should not be taken as a guide to the 
relative accuracy of the LB algorithm. 

3.2 Radial distribution functions (RDFs) 

In Figures [l] and [2] we plot the radial distribution function g{r)^ and the projections of 
the 'molecular' pair distribution function onto rotational invariants [50l[5l], as measured in 
various LB runs: 




(4) 
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Method 




U'^/NkBT 





0.10 


529 


643 


LB 


- 


0.08671 ± 0.00028 






529 


643 


BD 


- 


0.0870 ± 0.0004 






529 


22156a3 


MC 


- 


0.08757 ± 0.0001 
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0.10 


529 


643 


LB 


-2.929 ± 0.003 


0.2923 ± 0.0006 






529 


643 


BD 


-2.964 ± 0.002 


0.2935 ± 0.0006 






529 


22156a3 


MC 


-2.8830 ± 0.0008 


0.2850 ± 0.0002 
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0.10 


529 


643 


LB 


-11.811 ±0.002 


1.1692 ±0.0007 






529 


643 


BD 


-11.609 ±0.002 


1.1253 ±0.0007 






529 


22156a3 


MC 


-11.565 ±0.003 


1.1196 ±0.0006 


4 


0.20 


8239 


1283 


LB 


-3.966 ± 0.001 


0.5140 ± 0.0004 






8239 


1283 


BD 


-3.902 ± 0.001 


0.4970 ± 0.0006 






529 


11079a3 


MC 


-4.1895 ±0.0008 


0.4534 ± 0.0002 


8 


0.20 


8239 


1283 


LB 


-11.833 ±0.003 


1.233 ±0.001 






8239 


1283 


BD 


-11.646 ±0.003 


1.188 ±0.002 






529 


11079a3 


MC 


-11.677 ±0.003 


1.1925 ±0.0006 



Table 1: Energy equilibration data for LB, BD, and MC simulation runs. The quoted 
statistical errors are estimated on the basis of one standard deviation. A is the dipolar 
coupling constant and (j) = Nv{)/V is the volume fraction where is the number of colloids 
and = 47ra^/3 is the volume of one colloid. The system volumes are reported in lattice 
units for LB and BD runs, while the MC volumes are reported in units of the hard-core 
radius a (equal to 2.3 in lattice units). 

\ i<j I 

h22o{r) = ^^J^(E'^(--^.)[3(^.-^.r-l])- (7) 

Note that all the runs reported have = 0.10. We do not report equilibrium structural 
data for = 0.20 here because, even though the energy can be equilibrated for V = 64^, 
the modest number of particles combined with a very long autocorrelation time means that 
good statistics cannot be gathered for structural quantities even with millions of timesteps. 
For V = 128^ the statistics are better due to larger A^, but the maximum available runtime 

10^) is not long enough to guarantee equilibration as discussed above. 

Each of the above RDFs could equally be plotted in Fourier space (with g{r) then trans- 
forming into the static structure factor S{q))^ but the real space versions offer the more 
sensitive tests of equilibration. This is because, for the reasons already discussed, any errors 
are likely to occur in a localized range of r at close contact (r 2a). We find excellent 
agreement between LB, BD, and MC for g{r) at A = (data not shown). Figure [l] shows 
adequate agreement between all methods at A = 4 although BD and LB both show slight 
discrepancies from the MC data (which should be the most accurate) in the neighborhood 
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1.5 2 

r/2a 

(a) g{r) 



(b) /iiio(r) 




1.5 2 

r/2a 



(c) hiuir) 



(d) /l220(^) 



Figure 1: Radial Distribution functions (a) g{r)^ (b) /iiio(r), (c) /iii2(^), and (d) /i22o(^) for 
A = 4 and = 0.10: (black circles) BD; (red squares) LB; (green diamonds) MC. 



of the first peak for each correlator. (A slight discrepancy for LB is also detectable near the 
first minimum of /iiio(^)-) 

For A = 8 - Figure [2] - there is a clear discrepancy between LB and the other two 
methods, overestimating by 10% or so the height of the first peak in all the RDFs. This 
error is consistent in sign and magnitude with the energy discrepancies in Table [l} and with 
the specific type of noise discretization error reported in section [2] Given the other sources 
of error in LB [22j, we consider it acceptable. 

Overall, the observed behavior of g{r) and hi-^i^m{r) is in good accord with that established 
in earlier studies |32l [5l]. With A = 0, g{r) shows only short-range correlations (data not 
shown); with A 1, the primary peak in g{r) becomes very strongly pronounced due to 
the high degree of particle association to form, at these low densities, chain-like aggregates. 
^110 (^) helps distinguish parallel from anti-parallel correlations between the dipole moments, 
^112 (^) contains (minus) the dipolar potential, and /i22o(^) picks out 'nematic' orientational 
ordering. All of these functions show positive peaks at short range - at intervals close to 
2a - confirming the prevalence of the 'nose-to-tail' parallel conformations of nearby dipole 
moments within chains. As expected, the peaks become more pronounced with increasing A. 
For a visual confirmation of the chaining. Figure |3] shows two snapshots from equilibrated LB 
simulations of 529-particle systems with A = 4 and A = 8 at = 0.10. Each particle is color- 
coded to refiect the total number of particles in the cluster to which it belongs; monomers, 
clusters with n = 2-4 particles, and clusters with n > 5 particles are given unique colors. (See 
section [K2j) The chain-like structural motif is clearly visible in both systems. Figures |l]and 
[2] show that interparticle correlations are short-ranged (compared to the box dimensions), 
and hence the thermodynamic properties should not show any pronounced finite-size eflFects. 
Nonetheless, the cluster network in Figure [3] appears to span the simulation cell, and so 
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2 3 

r/2a 



(a) g{r) 



(c) hiuir) 



r/2a 

(b) /iiio(r) 




r/2a 

(d) /l220(^) 



Figure 2: Radial distribution functions (a) g{r)^ (b) /iiio(^), (c) /iii2(^), and (d) /i22o(^) for 
A = 8 and = 0.10: (black circles) BD; (red squares) LB; (green diamonds) MC. 



we might anticipate some finite-size effects in the long-time, long- wavelength dynamics. We 
have simulated the largest possible system sizes throughout. 




Figure 3: Snapshots from LB simulations of = 529 colloids at a volume fraction (f) = 0.10: 
(a) A = 4; (b) A = 8. Each particle is color-coded to refiect the total number of particles in 
the cluster to which it belongs: (dark blue) monomers; (light blue) dimers; (green) trimers; 
(yellow) tetramers; (red) clusters with 5 or more particles. 
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4 Dynamic correlators in equilibrium 



We now present results for the intermediate scattering function 5'(g,t), and its orientational 
analogs. As for the static structure, we restrict attention to the case = 0.10 where we can 
combine complete equilibration with adequate statistical averaging. We examined relaxations 
in the Fourier components of the number density and magnetization density at wavevectors q 
commensurate with the periodic boundary conditions; to improve the statistics, we averaged 
the appropriate correlators at wavevectors of the same magnitude q = \q\. 



4.1 Intermediate scattering function S{q^t) 

Figure [4] shows S{q^t) (eq[T]) as a function of t/rjj for A = 0, 4, and 8 at = 0.10. In each 
case curves are plotted for three different wavevectors; one close to the peak {q = q^) in 
the static structure factor S{q)^ one larger, and one smaller. These linear-linear plots allow 
one to see clearly the long-time relaxation of the structure. (The short time dynamics is 



considered in section 4.3 below.) 




Figure 4: Intermediate scattering functions with (a) A = 0, (b) A = 4, and (c) and (d) A = 8 
with linear and logarithmic abscissas, respectively: (solid lines) BD; (dashed lines) LB. In 
each case the black lines are for qa = 1.1514 and the green lines are for qa = 4.0456. The 
red lines are for a wavevector close to the peak of S{q)^ as follows: (a) qa = 2.6139 at A = 0; 
(b) qa = 3.0715 at A = 4; (c) and (d) qa = 3.2409 at A = 8. 



In all cases, the effect of hydrodynamic interactions is to slow down the relaxation of 
S{q^t)] this effect is most marked for wavenumbers well below (which in any case relax 
more slowly than those at the peak) . The long-time relaxations are not far from exponential 
in all cases, and in particular show no sign of decomposition into separate a and (3 relaxation 
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processes as expected in colloidal systems on approach to a glass transition [T5] . Absence of 
the latter is confirmed by plotting the time on a logarithmic scale, as presented for A = 8 
in Figure |4](d). It is notable, however, that the slowing by hydrodynamic interactions of the 
long-time relaxation, at least for qa 1, is much diminished at strong dipolar interactions 
(A = 8). This might be taken as evidence that structural rearrangement in this case is 
controlled mainly by the energetics of aggregate rearrangement (breaking and reformation of 
dipolar contacts) , and is no longer limited by the rate at which solvent can fiow around the 
evolving structure - a state of affairs generally expected to hold for glassy colloids. Caution 
is needed before drawing such a conclusion, however; many authors would, on adopting this 
reasoning, expect BD and LB curves to superpose only after rescaling of time by the short- 
time diflFusion constant at the peak, ^J. As discussed in section 4.3, hydrodynamic 
interactions continue to cause a factor of 2 change in this quantity even for A = 8. 



4.2 Orientational relaxation 

Defining a wavevector-dependent dipole density M{q^ t) = Z^jLi exp [—iq • rj{t)\ we can 
construct orientational correlators from the longitudinal (L) and transverse (T) components 
Mi^ = {M • q)q and Mt = M - Mi^ [53j: 

F{q,t) = N-\M{q,t)-M{-q,0)) (8) 
F^{q,t) = N-\M^{q,t) • M^{-q,0)) (9) 
F^{q,t) = N-\M^{q,t)'M^{-q,Q)). (10) 

Data for Fi^{q^ t) and F^{q^ t) dit q are plotted in Figure [5] for A = 0, 4, and 8. For clarity, 
we omit F(g, t) since this is a simple average of the longitudinal and transverse parts. The 
results are broadly comparable to the relaxation of S{q^t) at similar wavevectors. We note 
that the longitudinal relaxations are slower than the transverse ones. This might be ascribed 
to the slow rotational diffusion of chain orientations with respect to the wave vector qr, as 
compared to faster librational motions of dipoles perpendicular to the local chain orientation. 
Figure [5](d) shows (for A = 8) the q dependence of F{q^t). This is again comparable to that 
for the density relaxation. Note, though, that M is not a conserved quantity and therefore, 
unlike the density, is not compelled to relax slowly for qa < 1. The fact that it does so 
suggests that M is enslaved to slow particle rearrangements, as would arise if the dipole 
moments inside a cluster were to adopt frozen orientations relative to the positions of the 
constituent particles over the cluster's lifetime. With a dipolar bonding energy of 
for two linearly aligned dipoles, such behavior is quite plausible. 



4.3 Short time diffusion 

The shape of S{q^ t) is partly characterized by a short-time collective diflFusion constant 



Ds{q) = - 



1 



dhiS{q,t) 
Jt 



111 



where [..]s denotes a measurement taken at time scales long enough that a single particle 
indeed moves diffusively, but short enough that its average displacement remains small com- 
pared to a (or, if smaller, the surface-to-surface separation from neighboring particles). For 
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Figure 5: Orientational relaxations resolved in to longitudinal (Fl) and transverse (Ft) 
correlation functions, at = 0.10 and (a) A = 0, (b) A = 4, and (c) and (d) A = 8. In (a), 
(b), and (c), black lines are BD and green lines are LB: (solid lines) Fl; (dotted lines) Ft. 
In (d), solid lines are BD and dashed lines are LB: (black lines - upper) qa = L1514; (red 
lines - middle) qa = 3.2409; (green lines - lower) qa = 4.0456. 

an isolated particle we therefore require Tj^^t^ <^t <^ td where we recall that = a^p/vj is 
the time scale for steady fluid motion to be established at the particle scale, = m/Gnrja is 
the velocity autocorrelation time of the particle (of mass m), and td = cl^/Dq is the time for 
a particle to diflFuse its own radius. For particles with caging or bonding at surface-to-surface 
separations /i, the requirement t <C r/^) is replaced by t ^ ruih /ay. Deflning 

one finds that in the absence of hydrodynamic interactions the 'hydrodynamic factor' H{q) 
is always unity for all g, whereas experiments on, e.g., hard sphere colloids show values that 
are not only smaller but also g-dependent [HI [54] • For example, in hard-sphere colloids, 
0.2 < H{q) < 0.6 at ^ 0.3 and 1 < < 4 [Ml EH, while H{q') ^ 0.8 at ^ 0.10. 

The numerical evaluation of Ds{q)^ and hence of H{q)^ carries significant difficulties as- 
sociated with finite size corrections [36j. That is, the long range nature of the hydrodynamic 
interactions, in conjunction with periodic boundary conditions, makes Ds{q^ N) very slow to 
converge with system size or equivalently with particle number = (I)V/vq at fixed (f). For 
the case of hard spheres, at least, this can be brought under good control at a semi-empirical 
level by adopting the following correction [3f 



DM ^ Ds{q,N) (rj^ 
Do Do I 7] 



1.7601 



N N 



(13) 
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Figure 6: Short time decay of ln[S{q^t) / S{q)] as a function of t for (a) A = 0, (b) A = 4, 
and (c) A = 8, showing the extent of the hnear regime in each case. The sohd hues are BD 
and the dashed hues are LB. 



where t^oo is the so-caUed high-frequency viscosity of the suspension and rj is the solvent 
viscosity. 



To evaluate eq[T3]we numerically measured the high frequency viscosities for fully equili- 
brated systems with the given repulsive short-range potential and dipolar long-range interac- 
tion, using the recipe by Ladd [Ml- This amounts to calculating the integrated stress-stress 
correlation in a time window that is long enough to relax fluid degrees of freedom but too 
short for the colloids to move signiflcantly. For A = 0, 4, and 8 we found rj^Q / r] io be 1.0532, 
1.0717, and 1.1687, respectively. A similar procedure was used in jlSj for the case of colloids 



with long range coulombic repulsions. However, since eq [T3| was invented to account for the 
observed flnite-size behavior of systems of hard spheres, its use in other systems remains 
empirically questionable. Below we therefore present data both for Ds{q^N) as actually 



measured and for Ds{q) as estimated via eq 13, but use the latter value to calculate H{q). 

The above caveat applies particularly when long-range (e.g., dipolar) interactions are 
present. Arguably such interactions should create their own flnite size corrections, some- 
what akin to those from hydrodynamics. In this case one might expect that, even with 
hydrodynamics switched off, the measured H{q) would show size-dependent deviations from 
unity. In the data reported below we indeed flnd H{q) values signiflcantly less than unity 
for BD at large A; however, we know of no method to correct for this and make no attempt 
to do so. 

Figure [6] shows representative {q g*) short time S{q^t) data for the three values of 
A studied at = 0.10. In accordance with expectation, the regime of short time diflFusion 
is established beyond a few hundred timesteps, and for A = and 4 there is thereafter a 
wide region of exponential decay within which Ds{q) can be measured easily. For A = 8 this 
window is foreshortened ^ which is not surprising since the short time regime should end on 
the timescale of particle collisions. (For high interaction strengths, particles are bonded to 
neighbors with which they collide frequently.) Nonetheless a reasonable numerical estimate 
of the decay rate Ds{q^ N) can be made. In practice this was done by flrst identifying by eye 
the time window for short-time diffusion and then fltting to the log-linear plot within this 
window at each q. Finally the data for distinct q values were binned (each bin containing 
roughly ten wavevectors) and the statistical error then estimated for the binned data. 

Figures [tJ |8| and|9]show plots of S{q)^ DQ/Ds{q)^ and H{q) generated from our dynamic 
datasets for A = 0, 4, and 8, respectively. We also show, for comparison, the uncorrected 
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Figure 7: Structural and diffusion data for A = and (j) = 0.10: (black circles) BD; (red 
squares) LB; (green line) MC. In (c) the upper dataset (green) is with the uncorrected 
Ds{q^ N)^ and the lower dataset (red) is with the corrected Ds{q). 



Do/Ds{q^ N) curves; S{q) data generated from MC to check accuracy; and direct comparison 
with our BD results for Ds{q) and H{q). For the BD data no finite-size correction was 
made; for A = we recover H{q) = 1 to simulation accuracy, with smaller values at larger 
A presumably attributable to finite size effects in the thermodynamic sector, as discussed 
above. (It is possible that, were these to be corrected, the H{q) curves for LB could depend 
less strongly on A than in the results shown here.) 

When hydrodynamics is switched on, we obtain values of H{q^) 0.6 — 0.8 for all three 
A values that are comparable to previous data on hard sphere colloids at = 0.10 [TBI El]- 
However, for large A, S{q) and H{q) both suggest a rising trend at small wavenumbers 
{q < q'*/3). The rise in S{q) at low q is consistent with the formation of large dipolar 
clusters (and specifically, chains [31]). Clustering should also reduce the hydrodynamic 
friction per particle ^ as is familiar from the fact that large clusters sediment more quickly 
under gravity. That is, the body force increases linearly with particle number n whereas 
the viscous friction scales with hydrodynamic radius which is generally sublinear. This 
reduction is consistent with the observed rise in H{q) at low q. Note however that this rise, 
although detectable beyond the scatter in the H{q) data itself, is sensitive to the treatment 
of finite-size corrections and is therefore provisional. 
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Figure 8: Structural and diffusion data for A = 4 and = 0.10: (black circles) BD; (red 
squares) LB; (green line) MC. In (c) the upper dataset (green) is with the uncorrected 
Ds{q^ N)^ and the lower dataset (red) is with the corrected Ds{q). 
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Figure 9: Structural and diffusion data for A = 8 and = 0.10: (black circles) BD; (red 
squares) LB; (green line) MC. In (c) the upper dataset (green) is with the uncorrected 
Ds{q^ N)^ and the lower dataset (red) is with the corrected Ds{q). 
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5 Transient dynamics and cluster formation 



We now consider the evolution of the structure foUowing an initial quench, at time zero, 
from an equilibrium state with A = to one with either A = 4 or A = 8. We have studied 
such quenches at = 0.03,0.10, and 0.20. 



5.1 Transient dipolar energy 



Figure [10] shows the relaxation of the dipolar energy for each volume fraction as a function 
of time t following the quench; LB and BD data are directly compared. In all cases, the 
effect of hydrodynamic interactions is to slow the approach to equilibrium. However, the 
effect is quite modest, and comparable to that reported earlier for S{q^t) in equilibrium. 
The relaxation time is increased by no more than roughly a factor two, even for A = 8. 
Note that addition of hydrodynamics is by no means guaranteed to slow down, rather than 
speed up, the approach to equilibrium. A familiar counterexample is binary fluid phase 
separation, where fluid flow of the two species creates a less dissipative, and hence faster, 
phase-separation route than pure diffusion at intermediate and late times [20l EH EH] . 




H -4 




-I2h 



Figure 10: Relaxation of the dipolar energy following a quench from A = to (a) A = 4 and 
(b) A = 8: (black lines) BD; (orange lines) LB. Pairs of BD/LB curves correspond to the 
volume fractions, from top to bottom, cj) = 0.03, 0.10, and 0.20. 



5.2 Cluster statistics 

We define two dipolar particles to be in a bonded configuration if their pair dipolar interaction 
energy from eq [s] obeys 

[/^ < -0.75A. (14) 

This definition is somewhat arbitrary: in principle once could choose either an energy-based 
or a geometric criterion. Our choice corresponds to an energy criterion set by an equipotential 
surface, in configuration space, of the dipolar part of the interaction. This is more suitable 
than a criterion based solely on r: the latter would count as a bond any close encounter 
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Figure 11: Relaxation of cluster probabilities Pn{t) following quenches from A = to A = 4 
and 8 at = 0.03: (black circles) BD with A = 4; (green squares) LB with A = 4; (blue 
circles) BD with A = 8; (orange squares) LB with A = 8. 



between dipoles even if their orientation was such as to create a strongly repulsive force. In 
addition, our choice is designed to capture end-to-end bonding but reject most encounters 
between antiparallel dipoles even when their orientation is such as to create a bond. (Because 
of the short-range repulsion, the energy minimum for such bonds lies above the threshold in 
eq[l4j) The particular value of the energy threshold - which is intermediate between those 
used in earlier studies [55"^ ^ - gave cluster distributions in good accord with what was 
expected from visual inspection of simulation snapshots, and was sufficient for the current 
purpose of examining transient cluster formation. 



Using eq[T4]we partition each configuration of particles into a set of disjoint clusters, 
and monitor the fraction of particles that are assigned to clusters of size n. The time 
evolution of Pn{t) gives information about the growth of clusters following the quench from 
A = at t = 0. Figures 11 and 12 show Pn{t) data for various A at volume fractions = 0.03 
and 0.20. Once again, BD data is included for comparison. (The data is binned timewise 
with a stride of 25 timesteps for t < 200, 000 and 50 timesteps thereafter. This choice offered 
the best compromise between smoothness and sensitivity. The actual numbers of clusters 
are low, and the relative fiuctuations are high, so it is not easy to iron out the noise.) The 
transient Pn{t) dynamics is subject to a similar slowing by hydrodynamic interactions as was 
the energy transient. Other than this there are no obvious differences between the LB and 
BD data. For large A, both show characteristically peaked plots for P2, ^3, and P4 as small 
clusters build up and are then subsumed into larger ones. 
Finally, in Figure 
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we show the mean number of particles per cluster, A^^, as a function 
of time for the same runs. The same hydrodynamic slowing is evident. Mean cluster sizes 
in the LB simulations are slightly higher than those in the BD simulations. This can be 
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Relaxation of cluster probabilities Pn{t) following quenches from A = to A = 4 



0.20. Symbols as in Figure 11 



traced back to the small errors in dealing with particles close to contact, leading to more 
pronounced near-neighbour correlations, as discussed in section |3) 

6 Summary and Conclusions 

In this paper we have presented results for the equilibrium and transient dynamics of dipolar 
colloids with many-body hydrodynamic interactions. These were gained by incorporating the 
Ewald summation for the long-range dipolar interactions into our existing lattice Boltzmann 
algorithms, which handle the hydrodynamic forces by explicit propagation of momentum 
across the fluid residing on a lattice. The colloidal particles themselves are off-lattice and 
undergo molecular dynamics; Brownian motion is caused by fluctuating momentum transfer 
from the surrounding solvent, which creates correlated noise of the kind demanded by the 
fluctuation-dissipation theorem for hydrodynamically interacting particles. The full, fluid- 
driven noise can easily be replaced by local noise with no correlation between particles, 
creating a BD code. The resulting comparison of the LB and BD results allows the eflFect of 
many-body hydrodynamics to be isolated. 

At the volume fractions (0.03 < < 0.20) and interaction strengths (A = 4, 8) studied 
here, these eflFects are easily measurable but remain relatively modest. Quantitative shifts in 
both short-time and long-time diffusion were observed for wavevectors near and below the 
peak in the static structure factor. Likewise, we found shifts in transient relaxation rates for 
the cluster size distribution on approach to steady state following a quench from A = 0. In 
all cases, the system with hydrodynamic interactions relaxes more slowly than the equivalent 
BD system. However, for the range of volume fraction and interaction strengths studied. 
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Figure 13: Time evolution of mean cluster size for (a) A = 4 and (b) A = 8, at = 
0.03,0.10,0.20: (black circles) BD with = 0.03; (green squares) LB with = 0.03; (blue 
circles) BD with (j) = 0.10; (orange squares) LB with (j) = 0.10; (red circles) BD with = 0.20; 
(cyan squares) LB with = 0.20. 



this slowing down was rarely by more than a factor of two. 

Although it is possible that stronger hydrodynamic effects would be observed in the 
dynamics of a quiescent system at larger and A, their effects on long-time relaxation 
appear already to be decreasing for the largest values studied here. This could be a precursor 
to entering a glassy regime in which the crossing of local energy barriers limits relaxation 
rates; within this regime, conventional wisdom holds that hydrodynamics affects relaxational 
dynamics only through a scale factor [52]. However, despite the observation of slow transients 
and the difficulty of attaining full equilibration, even at = 0.20 and A = 8, we find no direct 
evidence for a glassy regime; specifically we see no separate a and /? relaxation processes. 
This is not surprising as our simulations run for at most 20 — 30r2:), and a truly glassy system 
would certainly not approach equilibrium, as ours do, on this time scale. 

Within our LB framework, there are serious obstacles to achieving much longer physical 
timescales using reasonable computational resources. One bottleneck remains the accurate 
treatment of noise; currently this requires a very large separation (order 10^ in our runs) 
between the simulation timestep and r^. Future algorithmic work will, we hope, partially 
address this issue. 

Our simulations were engineered to avoid the very large hydrodynamic forces that arise 
when hard colloidal particles come into lubrication contact. This was done by including 
a soft-core repulsion to maintain adequate separation between particle surfaces even when 
strongly bonded by the dipolar interactions. It is possible that these lubrication effects could 
enhance the relative role of hydrodynamic interactions, by further slowing the timescale 
for bond breakage and re-formation. To address this effect specifically (in a bulk periodic 
system) an algorithm such as ASD might be more suitable than LB. Note that within LB 
one can include a routine to address lubrication forces via an SD-like algorithm, but the 
computational scaling becomes bad when there are large clusters of particles in mutual 
lubrication contact. We do not know how well ASD would perform under such conditions, 
as compared to the purely repulsive interactions studied in [16]. 



19 



Even without lubrication, the effects of many-body hydrodynamics could, of course, also 
become much more pronounced in various nonequilibrium situations. These include the 
rheological response to steady and/or time-dependent shearing, and perhaps the nonlinear 
response to large orienting fields. We hope to address one or more of these topics in future 
work. 
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